%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Integrating GPE by the Method of Steepest Descent.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;

dt = .01;
g  =100;																	%Self-energy-ratio
R  =[-1:.05:1]*1.5;				 Z=R;
[r,z]=meshgrid(R,Z);
Psi=exp(-r.^2-z.^2);		    	Psi=Psi/norm(Psi);						%Trial-Wavefunction
%V  =r.^2+sqrt(8)*z.^2;														%Potential
V  =2.5*(z.^4-2*z.^2+1) + 5*r.^2;

for k=1:200
  [ gradPsi,Dumm]=gradient(Psi);
  [nablaPsi,Dumm]=gradient(gradPsi);
  H=-nablaPsi + V + g*Psi.^2;
  Psi=Psi.*(1-H*dt);
  Psi=Psi/norm(Psi);

figure(1);		clf;
subplot(221),surf(R,Z,V)
title('\fontsize{18}External Potential');
set(gca,'TickLength',[.03 .05],'FontSize',16,'LineWidth',1.2);
xlabel('\fontsize{18}{\itr/a_{ho}}');
ylabel('\fontsize{18}{\itz/a_{ho}}');
zlabel('\fontsize{18}{\itV/h\omega}');
subplot(222),surf(R,Z,Psi)
%axis([min(x(1,:)) max(x(1,:)) 0 5]);
title('\fontsize{18}BEC-Wavefunction');
set(gca,'TickLength',[.03 .05],'FontSize',16,'LineWidth',1.2);
xlabel('\fontsize{18}{\itr/a_{ho}}');
ylabel('\fontsize{18}{\itz/a_{ho}}');
zlabel('\fontsize{18}{\it\Psi/a_{ho}^{3/2}}');
pause(.01);
end;
